home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / cpoly.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  37.5 KB  |  1,045 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module cpoly)
  13.  
  14. ;;;This is a lisp version of algorithm 419 from the Communications of
  15. ;;;the ACM (p 97 vol 15 feb 1972) by Jenkins and Traub.
  16. ;;;That algorithm is followed very closely.
  17. ;;;Note the following modifications: arrays are indexed from 0 instead
  18. ;;;of 1.  This means that the variables n and nn are one less than the 
  19. ;;;acm verson.  The zeros are put into the arrays pr-sl and pi-sl, rather than
  20. ;;;into their own arrays.  The algorithm seems to benefit be taking are 
  21. ;;;mre 0.01 times the published values.
  22.  
  23. (declare-top
  24.  (*expr displa $listofvars meqhk displine)
  25.  (special logbas infin smalno are mre cr ci sr si tr ti zr zi n nn bool conv pvr
  26.       pvi $partswitch $keepfloat $demoivre $listconstvars $algebraic acp-sl
  27.       $polyfactor polysc polysc1 $ratfac $programmode)
  28.  (flonum logbas infin smalno are mre cr ci sr si tr ti zr zi xx yy cosr sinr bnd
  29.      xni t1 t2 otr oti svsr svsi pvr pvi mp ms omp relstp tp hvr hvi e ar
  30.      ai br bi x xm f dx df r1 $t hi lo max min acp-sl)
  31.  (fixnum degree n nn j l l1 l2 l3 cnt1 cnt2 jj polysc polysc1))
  32.  
  33. (declare-top(notype ($cpoly notype) (noshft-sl fixnum) (fxshft-sl fixnum)
  34.          (vrshft-sl fixnum) (calct-sl) (nexth-sl) (polyev-sl)
  35.          (cdivid-sl flonum flonum flonum flonum) (scale-sl))
  36.      (flonum (errev-sl flonum flonum) (cauchy-sl) (cmod-sl flonum flonum))
  37.      (fixnum (cpoly-sl fixnum))
  38.   #-PDP10(flonum (*f flonum flonum) (//f flonum flonum) (_f flonum fixnum))
  39.   #-PDP10(*expr *f //f _f)
  40.      (*lexpr $rat MAXIMA-ERROR)) 
  41.  
  42. ;;I changed these to unnamed arrays so as to be common lisp compatible
  43. ;;In older lisps use (defmacro aref (ar &rest dims) `(arraycall t ,@dims))
  44. ;;Really we should have non special variables for inner loop references to the
  45. ;;arrays, but.. later --wfs.
  46.  
  47. ;(declare 
  48. ; (array* (flonum *PR-SL* 1. *PI-SL* 1. *SHR-SL* 1. *SHI-SL* 1. *QPR-SL* 1. *QPI-SL* 1. *HR-SL* 1. *HI-SL*
  49. ;         1. *QHR-SL* 1. *QHI-SL* 1.)))
  50.  
  51. (declare-top
  52.  (special *PR-SL* *PI-SL* *SHR-SL* *SHI-SL* *QPR-SL* *QPI-SL* *HR-SL* *HI-SL*
  53.          *QHR-SL* *QHI-SL*))
  54. ;; Fixed for Twenex systems?
  55.  
  56. #+PDP10
  57. (and (not (get '*f 'subr)) 
  58.      (mapc '(lambda (x) (putprop x '(arith fasl dsk macsym) 'autoload))
  59.        '(*f //f +f -f _f)))
  60.  
  61. ;;; It is harder to underflow on lisp machine, but I suppose someday -BEE
  62.  
  63. #-PDP10
  64. (eval-when (compile load)
  65.        (defmacro *f  (a b) `(*$ ,a ,b))
  66.        (defmacro //f (a b) `(//$ ,a ,b))
  67.        (defmacro +f  (a b) `(+$ ,a ,b))
  68.        (defmacro -f  (a b) `(-$ ,a ,b))
  69.        )
  70.  
  71. ;(defmacro cpoly-large-flonum ()
  72. ;  #+LISPM (let ((a (float 0)))
  73. ;        (%p-dpb -1 1013 a) (%p-dpb -1 0007 a)
  74. ;        (%p-dpb-offset -1 0030 a 1) a)
  75. ;  #+NIL most-positive-double-float
  76. ;  #-(or LISPM NIL) '(fsc (lsh -1 -1) 0.)
  77. ;)
  78. ;
  79. ;(defmacro cpoly-small-flonum ()
  80. ;  #+LISPM (%p-dpb 100 0007 (float 0))
  81. ;  #+NIL least-positive-double-float
  82. ;  #-(or LISPM NIL) '(fsc (rot 1. -10.) 0.)
  83. ;  )
  84. ;
  85. ;(defmacro float-precision (pres)
  86. ;  pres ;Ignored on Lisp Machine
  87. ;  #+LISPM (let ((a (float 1)))
  88. ;        (%p-dpb-offset 1 0001 a 1) (f- a 1.0))
  89. ;  #+NIL (progn ;I wish i knew exactly what this meant. --gsb
  90. ;           `(*$ double-float-epsilon ,pres)
  91. ;           )
  92. ;  #-(or LISPM NIL) `(*$ (fsc (lsh 205. 26.) 0.) ,pres)
  93. ;  )
  94. ;;;
  95. #+Multics(defmacro _f  (a b) `(fsc ,a ,b))
  96.        
  97.  
  98. ;;; There were no comments on the following macros and functions, that on the
  99. ;;; Lisp Machine try to use %P-xxx subprimitives, and just used the #+LISPM
  100. ;;; conditional.  Of course, trying to transport this code to any other kind
  101. ;;; of Lisp Machine referneces random memory!  I've changed around the 
  102. ;;; conditionalization some, to make the macros and function only have definitions for
  103. ;;; machines that are known about, and I supported the 3600 as best as I could,
  104. ;;; since I had to intuit what the desired functionality is.  --HIC  3/5/83
  105.  
  106. ;;; Returns largest possible flonum?  (why not just a constant?)
  107. #-cl
  108. (defmacro cpoly-large-flonum ()
  109.       #-LISPM '(fsc (lsh -1 -1) 0.)
  110.       #+(and LISPM CADR)
  111.         (let ((a (float 0)))
  112.           (%p-dpb -1 1013 a) (%p-dpb -1 0007 a)
  113.           (%p-dpb-offset -1 0030 a 1) a)
  114.       #+(and LISPM 3600)
  115.         (times 2.0 1.45e38))  ;; was 2.9e38 but then file won't load 
  116.             ;; into non-Lispm systems.  E.g. on MC gives 
  117.             ;; NUMERIC-OVERFLOW - READ  - JPG
  118. ;; we are using long-float in maxima for flonum
  119. #+cl
  120. (defmacro cpoly-large-flonum () (float most-positive-long-float))
  121.  
  122. ;;; Returns smallest possible flonum?  (why not just a constant?)
  123. #-cl
  124. (defmacro cpoly-small-flonum ()
  125.       #-LISPM '(fsc (rot 1. -10.) 0.)
  126.       #+(and LISPM CADR) (%p-dpb 100 0007 (float 0))
  127.       #+(and LISPM 3600) 2.0e-38)
  128. #+cl
  129. (defmacro cpoly-small-flonum () (float least-positive-long-float))
  130.  
  131. #+clisp
  132. (eval-when (compile load eval)
  133.   (setq sys::*inhibit-floating-point-underflow* t))
  134.  
  135. ;;; Returns precision of floating point representation as least significant possible
  136. ;;; floating point number?  (again, why not just a constant?)
  137. #-cl
  138. (defmacro MAXIMA-float-precision (pres)
  139.       pres ;Ignored on Lisp Machine
  140.       #-LISPM `(*$ (fsc (lsh 205. 26.) 0.) ,pres)
  141.       #+(and LISPM CADR)
  142.         (let ((a (float 1)))
  143.           (%p-dpb-offset 1 0001 a 1)
  144.           (- a 1.0))
  145.       #+(and LISPM 3600)
  146.         (let* ((a (float 1))
  147.            (b (si:flonum-mantissa a nil))
  148.            (c (si:flonum-exponent a nil)))
  149.           (- (si:%flonum (dpb (f+ c (f+ 126. 24.)) si:%%float-exponent
  150.                   (dpb (logior b 1) si:%%float-fraction 0)))
  151.          1.0)))
  152. ;I doubt this is what was intended but I don't know --wfs.
  153.   #. (let ((e double-float-epsilon))
  154.        (cond ((= (float 1 e) (+ (float 1 e) e))
  155.           (format t
  156.               "Warning your double-float-epsilon is too small by CLtL"
  157.               ))))
  158.  
  159. #+cl
  160. (defun MAXIMA-float-precision (pres)
  161.   ;; Check double-float-epsilon is really ok for this lisp,
  162.   ;; according to the common lisp manual it should have.
  163.  
  164.      (cond ((eq (type-of pres) (type-of double-float-epsilon))
  165.         double-float-epsilon)
  166.        ((eq (type-of pres) (type-of short-float-epsilon))
  167.         short-float-epsilon)
  168.        (t (error "unknown type ~a" pres))))
  169.  
  170. ;;; _F scales its floating point argument by its scale argument, which is 
  171. ;;; a power of two to scale by.
  172.  
  173. #+cl
  174. (defun dfloat (x) (coerce x 'double-float))
  175.  
  176. ;; #+Franz
  177. ;; (defun _f (number scale) somebody-write-this)
  178.  
  179.  
  180. ;;; I don't think that this is called enough to justify making it quick.
  181. ;;; Ought to be faster--wfs
  182. #+(or Franz)
  183. (defun _f (number scale)
  184.     (cond ((zerop number) 0.0)
  185.       (t (f* number (cond ((zerop scale) 1)
  186.                  ((plusp scale)
  187.                   (^ 2.0 scale))
  188.                  ((// 1.0 (^ 2.0 (minus scale)))))))))
  189.  
  190. #+cl ;;should condition it to switch to double floats if too big
  191. (defun _f (number scale) (scale-float number scale))
  192.  
  193.  
  194. ;#+(and cl ti) ;;bug in microcoded
  195. ;(defun _f (f scale) (f* f (expt (float (float-radix f) f) scale)))
  196.  
  197.  
  198. ;#+(and LISPM CADR)
  199. ;(defun _f (number scale)            
  200. ;    (let ((ans (+ 0.0 number))            
  201. ;      (exp))
  202. ;      (setq exp (+ scale (%p-ldb 1013 ans)))    
  203. ;      (cond ((zerop number) 0.0)
  204. ;        ((> exp 3777) (error  "_F Overflow -- see MC:CFFK;CPOLY"))
  205. ;        ;; Should check zunderflow
  206. ;        ((< exp 0) (error  "_F Underflow -- see MC:CFFK;CPOLY"))
  207. ;        (t (%p-dpb exp 1013 ans) ans))))
  208.  
  209. ;#+(and LISPM 3600)
  210. ;(defun _f (number scale)
  211. ;  (let ((new-exp (f+ (f+ 126. 24.) (si:flonum-exponent number nil) scale))    
  212. ;    (mant (si:flonum-mantissa number nil)))
  213. ;    (cond ((zerop number) 0.0)
  214. ;      ((> new-exp 377) (error  "_F Overflow -- see MC:CFFK;CPOLY"))
  215. ;      ;; Should check zunderflow
  216. ;      ((< new-exp 0) (error  "_F Underflow -- see MC:CFFK;CPOLY"))
  217. ;      (t (si:%flonum (dpb new-exp si:%%float-exponent
  218. ;                  (dpb mant si:%%float-fraction 0)))))))
  219.  
  220. (setq acp-sl 0.2) 
  221.  
  222. (DEFMFUN $allroots (expr)
  223.        (prog (degree nn var res $partswitch $keepfloat $demoivre $listconstvars
  224.           $algebraic complex $ratfac den expr1)
  225.          (setq $keepfloat t $listconstvars t $algebraic t)
  226.          (setq expr1 (setq expr (meqhk expr)))
  227.          (setq var (delq '$%i (cdr ($listofvars expr))))
  228.          (or var (setq var (list (gensym))))
  229.          (cond ((not (= (length var) 1.))
  230.             (merror "ALLROOTS: polynomial not univariate: ~M" var))
  231.            ((setq var (car var))))
  232.          (setq expr ($rat expr '$%i var)
  233.            res (reverse (car (cdddar expr))))
  234.          (do ((i (f- (length res) (length (caddar expr))) (f1- i)))
  235.          ((= i 0.))
  236.          (setq res (cdr res)))
  237.          (setq den (cddr expr) expr (cadr expr))
  238.          ;;;check denominator is a complex number
  239.          (cond ((numberp den) (setq den (list den 0)))
  240.            ((eq (car den) (cadr res))
  241.             (setq den (cddr den))
  242.             (cond ((numberp (car den))
  243.                (cond ((null (cddr den)) (setq den (list 0 (car den))))
  244.                  ((numberp (caddr den))
  245.                   (setq den (list (caddr den) (car den))))
  246.                  (t (cpoly-err expr1))))
  247.               (t (cpoly-err expr1))))
  248.            (t (cpoly-err expr1)))
  249.          ;;;if the name variable has disappeared, this is caught here
  250.          (setq nn 0)
  251.          (cond ((numberp expr) (setq expr (list expr 0)))
  252.            ((eq (car expr) (car res)) (setq nn 1))
  253.            ((eq (car expr) (cadr res))
  254.             (setq expr (cddr expr))
  255.             (cond ((numberp (car expr))
  256.                (cond ((null (cddr expr)) (setq expr (list 0 (car expr))))
  257.                  ((numberp (caddr expr))
  258.                   (setq expr (list (caddr expr) (car expr))))
  259.                  (t (cpoly-err expr1))))
  260.               (t (cpoly-err expr1))))
  261.            (t (cpoly-err expr1)))
  262.          (cond ((= nn 0)
  263.             (cond ($polyfactor
  264.                ((lambda (cr ci)
  265.                     (cdivid-sl (float (car expr)) (float (cadr expr))
  266.                          (float (car den)) (float (cadr den)))
  267.                     (return (simplify (list '(mplus)
  268.                                 (simplify (list '(mtimes)
  269.                                         '$%i ci))
  270.                                 cr))))
  271.                 0.0 0.0))
  272.               (t (return (list '(mlist simp)))))))
  273.          (setq degree (cadr expr) nn (f1+ degree))
  274.          (setq *pr-sl* (*array nil 'flonum nn))
  275.          (setq *pi-sl* (*array nil 'flonum nn))
  276.          (or (catch 'notpoly
  277.             (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
  278.                     ((null expr))
  279.                     (setq l (f- degree (car expr)) res (cadr expr))
  280.                     (cond ((numberp res) (store (aref *PR-SL* l) (float res)))
  281.                       (t 
  282.                          (or (eq (car res) %i)
  283.                          (throw 'notpoly nil))
  284.                          (setq res (cddr res))
  285.                          (store (aref *PI-SL* l) (float (car res)))
  286.                          (setq res (caddr res))
  287.                          (and res (store (aref *PR-SL* l) (float res)))
  288.                          (setq complex t))))))
  289.          ;;;this should catch expressions like sin(x)-x
  290.          (progn (*rearray '*PR-SL*)
  291.             (*rearray '*PI-SL*)
  292.             (cpoly-err expr1)))
  293.          (setq *shr-sl* (*array nil 'flonum nn))
  294.          (setq *shi-sl* (*array nil 'flonum nn))
  295.          (setq *qpr-sl* (*array nil 'flonum nn))
  296.          (setq *hr-sl* (*array nil 'flonum degree))
  297.          (setq *qhr-sl*(*array nil 'flonum degree))
  298.          (cond
  299.            (complex
  300.          (setq *qpi-sl*(*array nil 'flonum nn))
  301.          (setq *hi-sl*(*array nil 'flonum degree))
  302.          (setq *qhi-sl*(*array nil 'flonum degree))
  303.                 ))
  304.          (setq nn degree)
  305.          (cond (complex (setq res (errset (cpoly-sl degree))))
  306.            ((setq res (errset (rpoly-sl degree)))))
  307.          (*rearray '*SHR-SL*)
  308.          (*rearray '*SHI-SL*)
  309.          (*rearray '*QPR-SL*)
  310.          (*rearray '*HR-SL*)
  311.          (*rearray '*QHR-SL*)
  312.          (cond (complex (*rearray '*QPI-SL*)
  313.                 (*rearray '*HI-SL*)
  314.                 (*rearray '*QHI-SL*)))
  315.          (or res
  316.          (mtell "~%Unexpected error. Treat results with caution."))
  317.          (cond ((= nn degree)
  318.             (*rearray '*PR-SL*)
  319.             (*rearray '*PI-SL*)
  320.             (merror "~%No roots found")))
  321.          (setq res nil)
  322.          (cond
  323.           ((not (zerop nn))
  324.            (mtell "~%Only ~S out of ~S roots found "
  325.               (f- degree nn) degree)
  326.            (setq expr 0.0)
  327.            (do
  328.         ((i 0. (f1+ i)))
  329.         ((> i nn))
  330.         (setq 
  331.          expr
  332.          (simplify
  333.           (list '(mplus)
  334.             expr
  335.             (simplify (list '(mtimes)
  336.                     (simplify (list '(mplus)
  337.                             (simplify (list '(mtimes)
  338.                                     '$%i
  339.                                     (aref *PI-SL* i)))
  340.                             (aref *PR-SL* i)))
  341.                     (simplify (list '(mexpt)
  342.                             var
  343.                             (f- nn i)))))))))
  344.            (setq res (cons expr res)))
  345.           ($polyfactor
  346.            (setq expr ((lambda (cr ci)
  347.                    (cdivid-sl (aref *PR-SL* 0) (aref *PI-SL* 0)
  348.                         (float (car den))
  349.                         (float (cadr den)))
  350.                    (simplify (list '(mplus)
  351.                            (simplify (list '(mtimes)
  352.                                    '$%i ci))
  353.                            cr)))
  354.                0.0 0.0)
  355.              res (cons expr res))))
  356.          (do
  357.           ((i degree (f1- i)))
  358.           ((= i nn))
  359.           (setq expr (simplify (list '(mplus)
  360.                      (simplify (list '(mtimes)
  361.                              '$%i
  362.                              (aref *PI-SL* i)))
  363.                      (aref *PR-SL* i))))
  364.           (setq 
  365.            res
  366.            (cond
  367.         ($polyfactor (cons (cond ((or complex (zerop (aref *PI-SL* i) ))
  368.                       (simplify (list '(mplus)
  369.                               var
  370.                               (simplify (list '(mminus)
  371.                                       expr)))))
  372.                      (t (setq i (f1- i))
  373.                         (simplify (list '(mplus)
  374.                                 (simplify (list '(mexpt)
  375.                                         var
  376.                                         2.))
  377.                                 (simplify (list '(mtimes)
  378.                                         var
  379.                                         (aref *PR-SL* i)))
  380.                                 (aref *PR-SL* (f1+ i))))))
  381.                    res))
  382.         ((cons ((lambda (expr) (cond ($programmode expr)
  383.                          (t (displine expr))))
  384.             (simplify (list '(mequal) var expr)))
  385.                res)))))
  386.          (*rearray '*PR-SL*)
  387.          (*rearray '*PI-SL*)
  388.          (return (simplify (cond ($polyfactor (cons '(mtimes) res))
  389.                      ((cons '(mlist) (nreverse res)))))))) 
  390.  
  391. (defun cpoly-err (expr)
  392.   (merror "ALLROOTS: not a polynomial:~%~M" expr))
  393.  
  394. (defun cpoly-sl (degree) 
  395.        ((lambda (logbas infin smalno are mre xx yy cosr sinr cr ci sr si tr ti zr zi bnd
  396.          n polysc polysc1 conv) 
  397.         (setq mre (*$ 2.0 (sqrt 2.0) are) yy (-$ xx))
  398.         (do ((i degree (f1- i)))
  399.             ((not (and (zerop (aref *PR-SL* i)) (zerop (aref *PI-SL* i)))) (setq nn i n (f1- i))))
  400.         (setq degree nn)
  401.         (do ((i 0. (f1+ i)))
  402.             ((> i nn))
  403.             (store (aref *SHR-SL* i) (cmod-sl (aref *PR-SL* i) (aref *PI-SL* i))))
  404.         (scale-sl )
  405.         (do nil
  406.             ((> 2. nn)
  407.              (cdivid-sl (-$ (aref *PR-SL* 1.)) (-$ (aref *PI-SL* 1.)) (aref *PR-SL* 0.) (aref *PI-SL* 0.))
  408.              (store (aref *PR-SL* 1.) cr)
  409.              (store (aref *PI-SL* 1.) ci)
  410.              (setq nn 0.))
  411.             (do ((i 0. (f1+ i)))
  412.             ((> i nn))
  413.             (store (aref *SHR-SL* i) (cmod-sl (aref *PR-SL* i) (aref *PI-SL* i))))
  414.             (setq bnd (cauchy-sl ))
  415.             (catch 'newroot
  416.                (do ((cnt1 1. (f1+ cnt1)))
  417.                    ((> cnt1 2.))
  418.                    (noshft-sl 5.)
  419.                    (do ((cnt2 1. (f1+ cnt2)))
  420.                    ((> cnt2 9.))
  421.                    (setq xx (prog2 nil
  422.                            (-$ (*$ cosr xx) (*$ sinr yy))
  423.                            (setq yy (+$ (*$ sinr xx)
  424.                                 (*$ cosr yy)))) 
  425.                      sr (*$ bnd xx) 
  426.                      si (*$ bnd yy))
  427.                    (fxshft-sl (f* 10. cnt2))
  428.                    (cond (conv (store (aref *PR-SL* nn) zr)
  429.                            (store (aref *PI-SL* nn) zi)
  430.                            (setq nn n n (f1- n))
  431.                            (do ((i 0. (f1+ i)))
  432.                            ((> i nn))
  433.                            (store (aref *PR-SL* i) (aref *QPR-SL* i))
  434.                            (store (aref *PI-SL* i) (aref *QPI-SL* i)))
  435.                            (throw 'newroot t))))))
  436.             (or conv (return t)))
  437.         (do ((i (f1+ nn) (f1+ i)))
  438.             ((> i degree))
  439.             (store (aref *PR-SL* i) (_f (aref *PR-SL* i) polysc1))
  440.             (store (aref *PI-SL* i) (_f (aref *PI-SL* i) polysc1)))
  441.         (do ((i 0. (f1+ i)) (j (f- polysc (f* polysc1 degree)) (f+ j polysc1)))
  442.             ((> i nn))
  443.             (store (aref *PR-SL* i) (_f (aref *PR-SL* i) j))
  444.             (store (aref *PI-SL* i) (_f (aref *PI-SL* i) j)))
  445.         nn)
  446.     (log 2.0) (cpoly-large-flonum)
  447.     (cpoly-small-flonum) (MAXIMA-float-precision acp-sl )
  448.     0.0 0.70710677 0.0 -0.069756474 0.99756405
  449.     0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0. 0. 0. nil)) 
  450.  
  451. (defun noshft-sl (l1) 
  452.        (do ((i 0. (f1+ i)) (xni (float nn) (1-$ xni)) (t1 (//$ (float nn))))
  453.        ((> i n))
  454.        (store (aref *HR-SL* i) (*$ (aref *PR-SL* i) xni t1))
  455.        (store (aref *HI-SL* i) (*$ (aref *PI-SL* i) xni t1)))
  456.        (do ((jj 1. (f1+ jj)))
  457.        ((> jj l1))
  458.        (cond ((> (cmod-sl (aref *HR-SL* n) (aref *HI-SL* n)) (*$ 10.0 are (cmod-sl (aref *PR-SL* n) (aref *PI-SL* n))))
  459.           (cdivid-sl (-$ (aref *PR-SL* nn)) (-$ (aref *PI-SL* nn)) (aref *HR-SL* n) (aref *HI-SL* n))
  460.           (setq tr cr ti ci)
  461.           (do ((j n (f1- j)) (t1) (t2))
  462.               ((> 1. j))
  463.               (setq t1 (aref *HR-SL* (f1- j)) t2 (aref *HI-SL* (f1- j)))
  464.               (store (aref *HR-SL* j) (-$ (+$ (aref *PR-SL* j) (*f t1 tr)) (*f t2 ti)))
  465.               (store (aref *HI-SL* j) (+$ (aref *PI-SL* j) (*f t1 ti) (*f t2 tr))))
  466.           (store (aref *HR-SL* 0.) (aref *PR-SL* 0.))
  467.           (store (aref *HI-SL* 0.) (aref *PI-SL* 0.)))
  468.          (t (do ((j n (f1- j)))
  469.             ((> 1. j))
  470.             (store (aref *HR-SL* j) (aref *HR-SL* (f1- j)))
  471.             (store (aref *HI-SL* j) (aref *HI-SL* (f1- j))))
  472.             (store (aref *HR-SL* 0.) 0.0)
  473.             (store (aref *HI-SL* 0.) 0.0))))) 
  474.  
  475. (defun fxshft-sl (l2) 
  476.        ((lambda (test pasd otr oti svsr svsi bool pvr pvi) 
  477.         (polyev-sl)
  478.         (setq conv nil)
  479.         (calct-sl)
  480.         (do ((j 1. (f1+ j)))
  481.             ((> j l2))
  482.             (setq otr tr oti ti)
  483.             (nexth-sl)
  484.             (calct-sl)
  485.             (setq zr (+$ sr tr) zi (+$ si ti))
  486.             (cond ((and (not bool) test (not (= j l2)))
  487.                (cond ((> (*$ 0.5 (cmod-sl zr zi))
  488.                      (cmod-sl (-$ tr otr) (-$ ti oti)))
  489.                   (cond (pasd (do ((i 0. (f1+ i)))
  490.                           ((> i n))
  491.                           (store (aref *SHR-SL* i) (aref *HR-SL* i))
  492.                           (store (aref *SHI-SL* i) (aref *HI-SL* i)))
  493.                           (setq svsr sr svsi si)
  494.                           (vrshft-sl 10.)
  495.                           (and conv (return nil))
  496.                           (setq test nil)
  497.                           (do ((i 0. (f1+ i)))
  498.                           ((> i n))
  499.                           (store (aref *HR-SL* i) (aref *SHR-SL* i))
  500.                           (store (aref *HI-SL* i) (aref *SHI-SL* i)))
  501.                           (setq sr svsr si svsi)
  502.                           (polyev-sl)
  503.                           (calct-sl))
  504.                     ((setq pasd t))))
  505.                  ((setq pasd nil))))))
  506.         (or conv (vrshft-sl 10.))
  507.         nil)
  508.     t nil 0.0 0.0 0.0 0.0 nil 0.0 0.0)) 
  509.  
  510. (defun vrshft-sl (l3) 
  511.        (setq conv nil sr zr si zi)
  512.        (do ((i 1. (f1+ i)) (bool1 nil) (mp) (ms) (omp) (relstp) (tp) (r1))
  513.        ((> i l3))
  514.        (polyev-sl)
  515.        (setq mp (cmod-sl pvr pvi) ms (cmod-sl sr si))
  516.        (cond ((> (*$ 20.0 (errev-sl ms mp)) mp)
  517.           (setq conv t zr sr zi si)
  518.           (return t)))
  519.        (cond ((= i 1.) (setq omp mp))
  520.          ((or bool1 (> omp mp) (not (< relstp 0.05)))
  521.           (cond ((> (*$ 0.1 mp) omp) (return t)) (t (setq omp mp))))
  522.          (t (setq tp relstp bool1 t)
  523.             (cond ((> are relstp) (setq tp are)))
  524.             (setq r1 (sqrt tp) 
  525.               sr (prog2 nil
  526.                     (-$ (*$ (1+$ r1) sr) (*f r1 si))
  527.                     (setq si (+$ (*$ (1+$ r1) si) (*f r1 sr)))))
  528.             (polyev-sl)
  529.             (do ((j 1. (f1+ j))) ((> j 5.)) (calct-sl) (nexth-sl))
  530.             (setq omp infin)))
  531.        (calct-sl)
  532.        (nexth-sl)
  533.        (calct-sl)
  534.        (or bool
  535.            (setq relstp (//$ (cmod-sl tr ti) (cmod-sl sr si)) 
  536.              sr (+$ sr tr) 
  537.              si (+$ si ti))))) 
  538.  
  539. (defun calct-sl nil 
  540.        (do ((i 1. (f1+ i))
  541.         ($t)
  542.         (hvr (store (aref *QHR-SL* 0.) (aref *HR-SL* 0.)))
  543.         (hvi (store (aref *QHI-SL* 0.) (aref *HI-SL* 0.))))
  544.        ((> i n)
  545.         (setq bool (not (> (cmod-sl hvr hvi) (*$ 10.0 are (cmod-sl (aref *HR-SL* n) (aref *HI-SL* n))))))
  546.         (cond ((not bool) (cdivid-sl (-$ pvr) (-$ pvi) hvr hvi) (setq tr cr ti ci))
  547.           (t (setq tr 0.0 ti 0.0)))
  548.         nil)
  549.        (setq $t (-$ (+$ (aref *HR-SL* i) (*f hvr sr)) (*f hvi si)))
  550.        (store (aref *QHI-SL* i) (setq hvi (+$ (aref *HI-SL* i) (*f hvr si) (*f hvi sr))))
  551.        (store (aref *QHR-SL* i) (setq hvr $t)))) 
  552.  
  553. (defun nexth-sl nil 
  554.        (cond (bool (do ((j 1. (f1+ j)))
  555.                ((> j n))
  556.                (store (aref *HR-SL* j) (aref *QHR-SL* (f1- j)))
  557.                (store (aref *HI-SL* j) (aref *QHI-SL* (f1- j))))
  558.            (store (aref *HR-SL* 0.) 0.0)
  559.            (store (aref *HI-SL* 0.) 0.0))
  560.          (t (do ((j 1. (f1+ j)) (t1) (t2))
  561.             ((> j n))
  562.             (setq t1 (aref *QHR-SL* (f1- j)) t2 (aref *QHI-SL* (f1- j)))
  563.             (store (aref *HR-SL* j) (-$ (+$ (aref *QPR-SL* j) (*f t1 tr)) (*f t2 ti)))
  564.             (store (aref *HI-SL* j) (+$ (aref *QPI-SL* j) (*f t1 ti) (*f t2 tr))))
  565.         (store (aref *HR-SL* 0.) (aref *QPR-SL* 0.))
  566.         (store (aref *HI-SL* 0.) (aref *QPI-SL* 0.))))
  567.        nil) 
  568.  
  569. (defun polyev-sl nil 
  570.        (setq pvr (store (aref *QPR-SL* 0.) (aref *PR-SL* 0.)) pvi (store (aref *QPI-SL* 0.) (aref *PI-SL* 0.)))
  571.        (do ((i 1. (f1+ i)) ($t))
  572.        ((> i nn))
  573.        (setq $t (-$ (+$ (aref *PR-SL* i) (*f pvr sr)) (*f pvi si)))
  574.        (store (aref *QPI-SL* i) (setq pvi (+$ (aref *PI-SL* i) (*f pvr si) (*f pvi sr))))
  575.        (store (aref *QPR-SL* i) (setq pvr $t)))) 
  576.  
  577. (defun errev-sl (ms mp) 
  578.        (-$ (*$ (do ((j 0. (f1+ j))
  579.             (e (//$ (*$ (cmod-sl (aref *QPR-SL* 0.) (aref *QPI-SL* 0.)) mre) (+$ are mre))))
  580.            ((> j nn) e)
  581.            (setq e (+$ (cmod-sl (aref *QPR-SL* j) (aref *QPI-SL* j)) (*$ e ms))))
  582.            (+$ are mre))
  583.        (*$ mp mre))) 
  584.  
  585. (defun cauchy-sl nil 
  586.        ((lambda (x xm) 
  587.         (store (aref *SHR-SL* nn) (-$ (aref *SHR-SL* nn)))
  588.         (cond ((not (zerop (aref *SHR-SL* n)))
  589.                (setq xm (-$ (//$ (aref *SHR-SL* nn) (aref *SHR-SL* n))))
  590.                (cond ((> x xm) (setq x xm)))))
  591.         (do ((f))
  592.             (nil)
  593.             (setq xm (*$ 0.1 x) f (aref *SHR-SL* 0.))
  594.             (do ((i 1. (f1+ i))) ((> i nn)) (setq f (+$ (aref *SHR-SL* i) (*f f xm))))
  595.             (cond ((not (< 0.0 f)) (return t)))
  596.             (setq x xm))
  597.         (do ((dx x) (df) (f))
  598.             ((> 5.0e-3 (abs (//$ dx x))) x)
  599.             (setq f (aref *SHR-SL* 0.) df f)
  600.             (do ((i 1. (f1+ i)))
  601.             ((> i n))
  602.             (setq f (+$ (*$ f x) (aref *SHR-SL* i)) df (+$ (*$ df x) f)))
  603.             (setq f (+$ (*$ f x) (aref *SHR-SL* nn)) dx (//$ f df) x (-$ x dx))))
  604.     (exp (//$ (-$ (log (aref *SHR-SL* nn)) (log (aref *SHR-SL* 0.))) (float nn)))
  605.     0.0)) 
  606.  
  607. (defun scale-sl nil
  608.        (do ((i 0. (f1+ i)) (j 0.) (x 0.0) (dx 0.0))
  609.        ((> i nn)
  610.         (setq x (//$ x (float (f- (f1+ nn) j)))
  611.           dx (//$ (-$ (log (aref *SHR-SL* nn)) (log (aref *SHR-SL* 0.))) (float nn))
  612.           polysc1 (fix (+$ 0.5 (//$ dx logbas)))
  613.           x (+$ x (*$ (float (f* polysc1 nn)) logbas 0.5))
  614.           polysc (fix (+$ 0.5 (//$ x logbas)))))
  615.        (cond ((zerop (aref *SHR-SL* i)) (setq j (f1+ j)))
  616.          (t (setq x (+$ x (log (aref *SHR-SL* i)))))))
  617.        (do ((i nn (f1- i)) (j (f- polysc) (f+ j polysc1)))
  618.        ((< i 0.))
  619.        (store (aref *PR-SL* i) (_f (aref *PR-SL* i) j))
  620.        (store (aref *PI-SL* i) (_f (aref *PI-SL* i) j))))
  621.  
  622. ;; (defun scale-sl nil 
  623. ;;        ((lambda (hi lo max min x l) 
  624. ;;         (do ((i 0. (f1+ i)))
  625. ;;             ((> i nn))
  626. ;;             (setq x (shr-sl i))
  627. ;;             (cond ((> x max) (setq max x)))
  628. ;;             (cond ((and (not (= x 0.0)) (< x min)) (setq min x))))
  629. ;;         (cond ((or (> lo min) (> max hi))
  630. ;;                (setq x (//$ lo min))
  631. ;;                (cond ((> x 1.0)
  632. ;;                   (cond ((> max (//$ infin x))
  633. ;;                      ;;;acm has < here but imsl agrees with me
  634. ;;                      (setq x 1.0))))
  635. ;;                  ((setq x (//$ (*$ (sqrt max) (sqrt min))))))
  636. ;;                (setq l (fix (+$ 0.5 (//$ (log x) logbas))))
  637. ;;                (cond ((not (= l 0.))
  638. ;;                   (do ((i 0. (f1+ i)))
  639. ;;                   ((> i nn))
  640. ;;                   (store (pr-sl i) (_f (pr-sl i) l))
  641. ;;                   (store (pi-sl i) (_f (pi-sl i) l)))))))
  642. ;;         l)
  643. ;;     (sqrt infin) (//$ smalno are) 0.0 infin 0.0 0.))
  644.  
  645. (defun cdivid-sl (ar ai br bi) 
  646.        ((lambda (r1) (cond ((and (zerop br) (zerop bi)) (setq cr (setq ci infin)))
  647.                ((> (abs bi) (abs br))
  648.                 (setq r1 (//f br bi) 
  649.                   bi (+$ bi (*f br r1)) 
  650.                   br (+$ ai (*f ar r1)) 
  651.                   cr (//f br bi) 
  652.                   br (-$ (*f ai r1) ar) 
  653.                   ci (//f br bi)))
  654.                ((setq r1 (//f bi br) 
  655.                   bi (+$ br (*f bi r1)) 
  656.                   br (+$ ar (*f ai r1)) 
  657.                   cr (//f br bi) 
  658.                   br (-$ ai (*f ar r1)) 
  659.                   ci (//f br bi)))))
  660.     0.0)
  661.        nil) 
  662.  
  663. (defun cmod-sl (ar ai) 
  664.        (setq ar (abs ar) ai (abs ai))
  665.        (cond ((> ai ar) (setq ar (//f ar ai)) (*$ ai (sqrt (1+$ (*f ar ar)))))
  666.          ((> ar ai) (setq ai (//f ai ar)) (*$ ar (sqrt (1+$ (*f ai ai)))))
  667.          ((*$ 1.41421357 ar)))) 
  668.  
  669. ;;*page
  670.  
  671. ;;;this is the algorithm for doing real polynomials.  it is algorithm 493 from
  672. ;;;acm toms vol 1 p 178 (1975) by jenkins.  note that array indexing starts from 0.
  673. ;;;the names of the arrays have been changed to be the same as for cpoly.
  674. ;;;the correspondence is:  p - pr-sl, qp - qpr-sl, k - hr-sl, qk - qhr-sl, svk - shr-sl,
  675. ;;;temp - shi-sl.  the roots are put in pr-sl and pi-sl.
  676. ;;;the variable si appears not to be used here
  677.  
  678. (declare-top(special sr u v a b c d a1 a3 a7 e f g h szr szi lzr lzi are mre n nn nz
  679.           #+gcl type ui vi s $polyfactor arp-sl)
  680.      (flonum a a0 a1 a3 a4 a5 a6 a7 aa are b b0 b1 b2 logbas bb betas betav bnd c c0
  681.          c1 c2 c3 c4 cc cosr d d0 e ee f g h infin kv lzi lzr mp mre ms omp
  682.          oss ots otv ovv pv relstp s sinr smalno sr ss svu svv szi szr t1 ts
  683.          tss tv tvv u ui v vi vv xx yy zm arp-sl)
  684.      (fixnum cnt degree i iflag j jj l l2 n nn nz #+gcl type)) 
  685.  
  686. (declare-top(fixnum (realit-sl))
  687.      (notype (rpoly-sl fixnum) (fxshfr-sl fixnum) (quadit-sl) (calcsc-sl) (nextk-sl)
  688.          (newest-sl) (quadsd-sl) (quad-sl flonum flonum flonum))) 
  689.  
  690. (setq arp-sl 1.0) 
  691.  
  692. (defun rpoly-sl (degree) 
  693.        ((lambda (logbas infin smalno are mre xx yy cosr sinr aa cc bb bnd sr u v t1 szr
  694.          szi lzr lzi nz n polysc polysc1 zerok conv1) 
  695.         (setq mre are yy (-$ xx))
  696.         (do ((i degree (f1- i))) ((not (zerop (aref *PR-SL* i))) (setq nn i n (f1- i))))
  697.         (setq degree nn)
  698.         (do ((i 0. (f1+ i))) ((> i nn)) (store (aref *SHR-SL* i) (abs (aref *PR-SL* i))))
  699.         (scale-sl)
  700.         (do nil
  701.             ((< nn 3.)
  702.              (cond ((= nn 2.)
  703.                 (quad-sl (aref *PR-SL* 0.) (aref *PR-SL* 1.) (aref *PR-SL* 2.))
  704.                 (cond ((and $polyfactor (not (zerop szi)))
  705.                    (store (aref *PR-SL* 2.) (//$ (aref *PR-SL* 2.) (aref *PR-SL* 0.)))
  706.                    (store (aref *PR-SL* 1.) (//$ (aref *PR-SL* 1.) (aref *PR-SL* 0.)))
  707.                    (store (aref *PI-SL* 2.) 1.0))
  708.                   (t (store (aref *PR-SL* 2.) szr)
  709.                      (store (aref *PI-SL* 2.) szi)
  710.                      (store (aref *PR-SL* 1.) lzr)
  711.                      (store (aref *PI-SL* 1.) lzi))))
  712.                (t (store (aref *PR-SL* 1.) (-$ (//$ (aref *PR-SL* 1.) (aref *PR-SL* 0.))))))
  713.              (setq nn 0.))
  714.             (do ((i 0. (f1+ i))) ((> i nn)) (store (aref *SHR-SL* i) (abs (aref *PR-SL* i))))
  715.             (setq bnd (cauchy-sl))
  716.             (do ((i 1. (f1+ i)))
  717.             ((> i n))
  718.             (store (aref *HR-SL* i) (//$ (*$ (float (f- n i)) (aref *PR-SL* i)) (float n))))
  719.             (store (aref *HR-SL* 0.) (aref *PR-SL* 0.))
  720.             (setq aa (aref *PR-SL* nn) bb (aref *PR-SL* n) zerok (zerop (aref *HR-SL* n)))
  721.             (do ((jj 1. (f1+ jj)))
  722.             ((> jj 5.))
  723.             (setq cc (aref *HR-SL* n))
  724.             (cond (zerok (do ((j n (f1- j)))
  725.                      ((< j 1.))
  726.                      (store (aref *HR-SL* j) (aref *HR-SL* (f1- j))))
  727.                      (store (aref *HR-SL* 0.) 0.0)
  728.                      (setq zerok (zerop (aref *HR-SL* n))))
  729.                   (t (setq t1 (-$ (//$ aa cc)))
  730.                  (do ((j n (f1- j)))
  731.                      ((< j 1.))
  732.                      (store (aref *HR-SL* j) (+$ (*$ t1 (aref *HR-SL* (f1- j))) (aref *PR-SL* j))))
  733.                  (store (aref *HR-SL* 0.) (aref *PR-SL* 0.))
  734.                  (setq zerok (not (> (abs (aref *HR-SL* n))
  735.                              (*$ (abs bb) are 10.0)))))))
  736.             (do ((i 0. (f1+ i))) ((> i n)) (store (aref *SHI-SL* i) (aref *HR-SL* i)))
  737.             (do ((cnt 1. (f1+ cnt)))
  738.             ((> cnt 20.) (setq conv1 nil))
  739.             (setq xx (prog2 nil
  740.                     (-$ (*$ cosr xx) (*$ sinr yy))
  741.                     (setq yy (+$ (*$ sinr xx) (*$ cosr yy)))) 
  742.                   sr (*$ bnd xx) 
  743.                   u (*$ -2.0 sr) 
  744.                   v bnd)
  745.             (fxshfr-sl (f* 20. cnt))
  746.             (cond ((> nz 0.)
  747.                    (store (aref *PR-SL* nn) szr)
  748.                    (store (aref *PI-SL* nn) szi)
  749.                    (cond ((= nz 2.)
  750.                       (store (aref *PR-SL* n) lzr)
  751.                       (store (aref *PI-SL* n) lzi)
  752.                       (cond ((and $polyfactor (not (zerop szi)))
  753.                          (store (aref *PR-SL* nn) v)
  754.                          (store (aref *PR-SL* n) u)
  755.                          (store (aref *PI-SL* nn) 1.0)))))
  756.                    (setq nn (f- nn nz) n (f1- nn))
  757.                    (do ((i 0. (f1+ i))) ((> i nn)) (store (aref *PR-SL* i) (aref *QPR-SL* i)))
  758.                    (return nil)))
  759.             (do ((i 0. (f1+ i))) ((> i n)) (store (aref *HR-SL* i) (aref *SHI-SL* i))))
  760.             (or conv1 (return nil)))
  761.         (cond ($polyfactor
  762.                (do ((i degree (f1- i)))
  763.                ((= i nn))
  764.                (cond ((zerop (aref *PI-SL* i))
  765.                   (store (aref *PR-SL* i) (_f (aref *PR-SL* i) polysc1)))
  766.                  (t (store (aref *PR-SL* i) (_f (aref *PR-SL* i) (f* 2. polysc1)))
  767.                     (setq i (f1- i))
  768.                     (store (aref *PR-SL* i) (_f (aref *PR-SL* i) polysc1))))))
  769.               (t (do ((i (f1+ nn) (f1+ i)))
  770.                  ((> i degree))
  771.                  (store (aref *PR-SL* i) (_f (aref *PR-SL* i) polysc1))
  772.                  (store (aref *PI-SL* i) (_f (aref *PI-SL* i) polysc1)))))
  773.         (do ((i 0. (f1+ i)) (j (f- polysc (f* polysc1 degree)) (f+ j polysc1)))
  774.             ((> i nn))
  775.             (store (aref *PR-SL* i) (_f (aref *PR-SL* i) j))))
  776.     (log 2.0) (cpoly-large-flonum)
  777.     (cpoly-small-flonum) (MAXIMA-float-precision arp-sl)
  778.     0.0 0.70710677 0.0 -0.069756474 0.99756405
  779.     0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0. 0. 0. 0. 0. t))  
  780.  
  781. (defun fxshfr-sl (l2) 
  782.        ((lambda (type a b c d e f g h a1 a3 a7) 
  783.      (setq nz 0.)
  784.      (quadsd-sl )
  785.      (calcsc-sl )
  786.      (do ((j 1. (f1+ j)) (betav 0.25) (betas 0.25)
  787.           (oss sr) (ovv v) (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
  788.           (ui) (vi) (s) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
  789.          ((> j l2))
  790.          (nextk-sl )
  791.          (calcsc-sl )
  792.          (newest-sl )
  793.          (setq vv vi ss 0.0)
  794.          (or (zerop (aref *HR-SL* n)) (setq ss (-$ (//$ (aref *PR-SL* nn) (aref *HR-SL* n)))))
  795.          (setq tv 1.0 ts 1.0)
  796.          (cond ((not (or (= j 1.) (= type 3.)))
  797.             (or (zerop vv) (setq tv (abs (//$ (-$ vv ovv) vv))))
  798.             (or (zerop ss ) (setq ts (abs (//$ (-$ ss oss) ss))))
  799.             (setq tvv 1.0)
  800.             (and (< tv otv) (setq tvv (*$ tv otv)))
  801.             (setq tss 1.0)
  802.             (and (< ts ots) (setq tss (*$ ts ots)))
  803.             (setq vpass (< tvv betav) spass (< tss betas))
  804.             (cond ((or spass vpass)
  805.                (setq svu u svv v)
  806.                (do ((i 0. (f1+ i))) ((> i n)) (store (aref *SHR-SL* i) (aref *HR-SL* i)))
  807.                (setq s ss vtry nil stry nil)
  808.                (and (do ((bool (not (and spass
  809.                              (or (not vpass) (< tss tvv))))
  810.                        t)
  811.                      (l50 nil nil))
  812.                     (nil)
  813.                     (cond (bool (quadit-sl )
  814.                         (and (> nz 0.) (return t))
  815.                         (setq vtry t betav (*$ 0.25 betav))
  816.                         (cond ((or stry (not spass))
  817.                                (setq l50 t))
  818.                               (t (do ((i 0. (f1+ i)))
  819.                                  ((> i n))
  820.                                  (store (aref *HR-SL* i)
  821.                                     (aref *SHR-SL* i)))))))
  822.                     (cond ((not l50)
  823.                        (setq iflag (realit-sl ))
  824.                        (and (> nz 0.) (return t))
  825.                        (setq stry t betas (*$ 0.25 betas))
  826.                        (cond ((zerop iflag) (setq l50 t))
  827.                          (t (setq ui (-$ (+$ s s)) 
  828.                               vi (*$ s s))))))
  829.                     (cond (l50 (setq u svu v svv)
  830.                            (do ((i 0. (f1+ i)))
  831.                            ((> i n))
  832.                            (store (aref *HR-SL* i) (aref *SHR-SL* i)))
  833.                            (and (or (not vpass) vtry)
  834.                             (return nil)))))
  835.                 (return nil))
  836.                (quadsd-sl )
  837.                (calcsc-sl )))))
  838.          (setq ovv vv oss ss otv tv ots ts)))
  839.     0. 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0)) 
  840.  
  841. (defun quadit-sl nil 
  842.        (setq nz 0. u ui v vi)
  843.        (do ((tried) (j 0.) (ee) (zm) (t1) (mp) (relstp) (omp))
  844.        (nil)
  845.        (quad-sl 1.0 u v)
  846.        (and (> (abs (-$ (abs szr) (abs lzr))) (*$ 0.01 (abs lzr))) (return nil))
  847.        (quadsd-sl )
  848.        (setq mp (+$ (abs (-$ a (*$ szr b))) (abs (*$ szi b))) 
  849.          zm (sqrt (abs v)) 
  850.          ee (*$ 2.0 (abs (aref *QPR-SL* 0.))) 
  851.          t1 (-$ (*$ szr b)))
  852.        (do ((i 1. (f1+ n))) ((> i n)) (setq ee (+$ (*$ ee zm) (abs (aref *QPR-SL* i)))))
  853.        (setq ee (+$ (*$ ee zm) (abs (+$ a t1))) 
  854.          ee (-$ (*$ (+$ (*$ 5.0 mre) (*$ 4.0 are)) ee)
  855.             (*$ (+$ (*$ 5.0 mre) (*$ 2.0 are))
  856.                 (+$ (abs (+$ a t1)) (*$ (abs b) zm)))
  857.             (*$ -2.0 are (abs t1))))
  858.        (cond ((not (> mp (*$ 20.0 ee))) (setq nz 2.) (return nil)))
  859.        (setq j (f1+ j))
  860.        (and (> j 20.) (return nil))
  861.        (cond ((not (or (< j 2.) (> relstp 0.01) (< mp omp) tried))
  862.           (and (< relstp are) (setq relstp are))
  863.           (setq relstp (sqrt relstp) 
  864.             u (-$ u (*$ u relstp)) 
  865.             v (+$ v (*$ v relstp)))
  866.           (quadsd-sl )
  867.           (do ((i 1. (f1+ i))) ((> i 5.)) (calcsc-sl ) (nextk-sl ))
  868.           (setq tried t j 0.)))
  869.        (setq omp mp)
  870.        (calcsc-sl )
  871.        (nextk-sl )
  872.        (calcsc-sl )
  873.        (newest-sl )
  874.        (and (zerop vi) (return nil))
  875.        (setq relstp (abs (//$ (-$ vi v) vi)) u ui v vi))) 
  876.  
  877. (defun realit-sl nil 
  878.        (setq nz 0.)
  879.        (do ((j 0.) (pv) (ee) (ms) (mp) (kv) (t1) (omp))
  880.        (nil)
  881.        (setq pv (aref *PR-SL* 0.))
  882.        (store (aref *QPR-SL* 0.) pv)
  883.        (do ((i 1. (f1+ i)))
  884.            ((> i nn))
  885.            (setq pv (+$ (*$ pv s) (aref *PR-SL* i)))
  886.            (store (aref *QPR-SL* i) pv))
  887.        (setq mp (abs pv) ms (abs s) ee (*$ (//$ mre (+$ are mre)) (abs (aref *QPR-SL* 0.))))
  888.        (do ((i 1. (f1+ i))) ((> i nn)) (setq ee (+$ (*$ ee ms) (abs (aref *QPR-SL* i)))))
  889.        (cond ((not (> mp (*$ 20.0 (-$ (*$ (+$ are mre) ee) (*$ mre mp)))))
  890.           (setq nz 1. szr s szi 0.0)
  891.           (return 0.)))
  892.        (setq j (f1+ j))
  893.        (and (> j 10.) (return 0.))
  894.        (cond ((not (or (< j 2.)
  895.                (> (abs t1) (*$ 1.0e-3 (abs (-$ s t1))))
  896.                (not (> mp omp))))
  897.           (return 1.)))
  898.        (setq omp mp kv (aref *HR-SL* 0.))
  899.        (store (aref *QHR-SL* 0.) kv)
  900.        (do ((i 1. (f1+ i)))
  901.            ((> i n))
  902.            (setq kv (+$ (*$ kv s) (aref *HR-SL* i)))
  903.            (store (aref *QHR-SL* i) kv))
  904.        (cond ((> (abs kv) (*$ (abs (aref *HR-SL* n)) 10.0 are))
  905.           (setq t1 (-$ (//$ pv kv)))
  906.           (store (aref *HR-SL* 0.) (aref *QPR-SL* 0.))
  907.           (do ((i 1. (f1+ i)))
  908.               ((> i n))
  909.               (store (aref *HR-SL* i) (+$ (*$ t1 (aref *QHR-SL* (f1- i))) (aref *QPR-SL* i)))))
  910.          (t (store (aref *HR-SL* 0.) 0.0)
  911.             (do ((i 1. (f1+ i))) ((> i n)) (store (aref *HR-SL* i) (aref *QHR-SL* (f1- i))))))
  912.        (setq kv (aref *HR-SL* 0.))
  913.        (do ((i 1. (f1+ i))) ((> i n)) (setq kv (+$ (*$ kv s) (aref *HR-SL* i))))
  914.        (setq t1 0.0)
  915.        (and (> (abs kv) (*$ (abs (aref *HR-SL* n)) 10.0 are)) (setq t1 (-$ (//$ pv kv))))
  916.        (setq s (+$ s t1)))) 
  917.  
  918. (defun calcsc-sl nil 
  919.        (setq d (aref *HR-SL* 0.))
  920.        (store (aref *QHR-SL* 0.) d)
  921.        (setq c (-$ (aref *HR-SL* 1.) (*$ u d)))
  922.        (store (aref *QHR-SL* 1.) c)
  923.        (do ((i 2. (f1+ i)) (c0))
  924.        ((> i n))
  925.        (setq c0 (-$ (aref *HR-SL* i) (*$ u c) (*$ v d)))
  926.        (store (aref *QHR-SL* i) c0)
  927.        (setq d c c c0))
  928.        (cond ((not (or (> (abs c) (*$ (abs (aref *HR-SL* n)) 100.0 are))
  929.                (> (abs d) (*$ (abs (aref *HR-SL* (f1- n))) 100.0 are))))
  930.           (setq type 3.))
  931.          ((not (< (abs d) (abs c)))
  932.           (setq type 2. 
  933.             e (//$ a d) 
  934.             f (//$ c d) 
  935.             g (*$ u b) 
  936.             h (*$ v b) 
  937.             a3 (+$ (*$ (+$ a g) e) (*$ h (//$ b d))) 
  938.             a1 (-$ (*$ b f) a) 
  939.             a7 (+$ (*$ (+$ f u) a) h)))
  940.          (t (setq type 1. 
  941.               e (//$ a c) 
  942.               f (//$ d c) 
  943.               g (*$ u e) 
  944.               h (*$ v b) 
  945.               a3 (+$ (*$ a e) (*$ (+$ (//$ h c) g) b)) 
  946.               a1 (-$ b (*$ a (//$ d c))) 
  947.               a7 (+$ a (*$ g d) (*$ h f)))))
  948.        nil) 
  949.  
  950. (defun nextk-sl nil 
  951.        (cond ((= type 3.)
  952.           (store (aref *HR-SL* 0.) 0.0)
  953.           (store (aref *HR-SL* 1.) 0.0)
  954.           (do ((i 2. (f1+ i))) ((> i n)) (store (aref *HR-SL* i) (aref *QHR-SL* (f- i 2.)))))
  955.          ((> (abs a1) (*$ (abs (cond ((= type 1.) b) (a))) 10.0 are))
  956.           (setq a7 (//$ a7 a1) a3 (//$ a3 a1))
  957.           (store (aref *HR-SL* 0.) (aref *QPR-SL* 0.))
  958.           (store (aref *HR-SL* 1.) (-$ (aref *QPR-SL* 1.) (*$ a7 (aref *QPR-SL* 0.))))
  959.           (do ((i 2. (f1+ i)))
  960.           ((> i n))
  961.           (store (aref *HR-SL* i)
  962.              (+$ (*$ a3 (aref *QHR-SL* (f- i 2.)))
  963.                  (-$ (*$ a7 (aref *QPR-SL* (f1- i))))
  964.                  (aref *QPR-SL* i)))))
  965.          (t (store (aref *HR-SL* 0.) 0.0)
  966.         (store (aref *HR-SL* 1.) (-$ (*$ a7 (aref *QPR-SL* 0.))))
  967.         (do ((i 2. (f1+ i)))
  968.             ((> i n))
  969.             (store (aref *HR-SL* i)
  970.                (-$ (*$ a3 (aref *QHR-SL* (f- i 2.))) (*$ a7 (aref *QPR-SL* (f1- i))))))))
  971.        nil) 
  972.  
  973. (defun newest-sl nil 
  974.        ((lambda (a4 a5 b1 b2 c1 c2 c3 c4) 
  975.         (cond ((= type 3.) (setq ui 0.0 vi 0.0))
  976.               (t (cond ((= type 2.)
  977.                 (setq a4 (+$ (*$ (+$ a g) f) h) 
  978.                       a5 (+$ (*$ (+$ f u) c) (*$ v d))))
  979.                    (t (setq a4 (+$ a (*$ u b) (*$ h f)) 
  980.                     a5 (+$ c (*$ (+$ u (*$ v f)) d)))))
  981.              (setq b1 (-$ (//$ (aref *HR-SL* n) (aref *PR-SL* nn))) 
  982.                    b2 (-$ (//$ (+$ (aref *HR-SL* (f1- n)) (*$ b1 (aref *PR-SL* n))) (aref *PR-SL* nn))) 
  983.                    c1 (*$ v b2 a1) 
  984.                    c2 (*$ b1 a7) 
  985.                    c3 (*$ b1 b1 a3) 
  986.                    c4 (-$ c1 c2 c3) 
  987.                    c1 (+$ a5 (*$ b1 a4) (-$ c4)))
  988.              (cond ((zerop c1) (setq ui 0.0 vi 0.0))
  989.                    (t (setq ui (-$ u
  990.                            (//$ (+$ (*$ u (+$ c3 c2))
  991.                             (*$ v
  992.                                 (+$ (*$ b1 a1)
  993.                                 (*$ b2 a7))))
  994.                             c1)) 
  995.                     vi (*$ v (1+$ (//$ c4 c1))))))))
  996.         nil)
  997.     0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0)) 
  998.  
  999. (defun quadsd-sl nil 
  1000.        (setq b (aref *PR-SL* 0.))
  1001.        (store (aref *QPR-SL* 0.) b)
  1002.        (setq a (-$ (aref *PR-SL* 1.) (*$ u b)))
  1003.        (store (aref *QPR-SL* 1.) a)
  1004.        (do ((i 2. (f1+ i)) (c0))
  1005.        ((> i nn))
  1006.        (setq c0 (-$ (aref *PR-SL* i) (*$ u a) (*$ v b)))
  1007.        (store (aref *QPR-SL* i) c0)
  1008.        (setq b a a c0))) 
  1009.  
  1010. (defun quad-sl (a0 b1 c0) 
  1011.        (setq szr 0.0 szi 0.0 lzr 0.0 lzi 0.0)
  1012.        ((lambda (b0 d0 e) 
  1013.         (cond ((zerop a0 ) (or (zerop b1 ) (setq szr (-$ (//$ c0 b1)))))
  1014.               ((zerop c0 ) (setq lzr (-$ (//$ b1 a0))))
  1015.               (t (setq b0 (//$ b1 2.0))
  1016.              (cond ((< (abs b0) (abs c0))
  1017.                 (setq e a0)
  1018.                 (and (< c0 0.0) (setq e (-$ a0)))
  1019.                 (setq e (-$ (*$ b0 (//$ b0 (abs c0))) e) 
  1020.                       d0 (*$ (sqrt (abs e)) (sqrt (abs c0)))))
  1021.                    (t (setq e (-$ 1.0 (*$ (//$ a0 b0) (//$ c0 b0))) 
  1022.                     d0 (*$ (sqrt (abs e)) (abs b0)))))
  1023.              (cond ((< e 0.0)
  1024.                 (setq szr (-$ (//$ b0 a0)) 
  1025.                       lzr szr 
  1026.                       szi (abs (//$ d0 a0)) 
  1027.                       lzi (-$ szi)))
  1028.                    (t (or (< b0 0.0) (setq d0 (-$ d0)))
  1029.                   (setq lzr (//$ (-$ d0 b0) a0))
  1030.                   (or (zerop lzr ) (setq szr (//$ (//$ c0 lzr) a0)))))))
  1031.         nil)
  1032.     0.0 0.0 0.0)) 
  1033.  
  1034. #-NIL
  1035. (declare-top(unspecial logbas infin smalno are mre cr ci sr si tr ti zr zi
  1036.             n nn bool conv pvr pvi acp-sl polysc polysc1 sr u v a
  1037.             b c d a1 a3 a7 e f g h szr szi lzr lzi are mre n nn nz
  1038.             #+gcl type ui vi s arp-sl ))
  1039.  
  1040. (declare-top(unspecial logbas infin smalno are mre cr ci sr si tr ti zr zi
  1041.             n nn bool conv pvr pvi acp-sl polysc polysc1 sr u v a
  1042.             b c d a1 a3 a7 e f g h szr szi lzr lzi are mre n nn nz
  1043.             #+gcl type ui vi s arp-sl))
  1044.  
  1045.